import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import py_tool
import matplotlib.dates as mdates
from datetime import datetime, timedelta
from matplotlib.ticker import AutoMinorLocator, NullFormatter, FixedLocator, FormatStrFormatter
import matplotlib.font_manager as font_manager
#从高震级开始标注
#先标注10 - 9, 9 - 8, 8 - 7
#尽量标注新地震;
#较近的地震只标注
def calpaintEqLabelIndex(SDATE, EDATE, EQDATE, EQRANK, EQPOS, EQDIS):
    N      = EDATE - SDATE
    Y      = int(N/365)
    Y3     = Y * 3 #有多少个地震都只绘制这么多地震;
    M      = len(EQRANK)
    d      = 200
    ind    = np.full(M, 0)    #如果序列为1,表示要绘制, 否者不绘制
    sorted_indices = np.argsort(-EQRANK)
    EQRANK         = EQRANK[sorted_indices]
    EQDATE         = EQDATE[sorted_indices]
    EQPOS          = EQPOS[sorted_indices]
    EQDIS          = EQDIS[sorted_indices]

    for i in range(len(ind)):
        if isANNO(EQDATE, i, d, ind) == 1 and i <= Y3:
            ind[i] = 1
    return EQRANK, EQDATE, EQPOS, EQDIS, ind

#判断列表中是否有比该地震更大的地震
def isANNO(EQDATE, selEQIndex, d, ind):
    status = 0 # 不画
    iEQDATE = EQDATE[selEQIndex]
    isExist = ExistNearEq(EQDATE, selEQIndex, d); # isExist = 0 附近没有地震; isExist = 1 附近有地震;
    for i in range(len(EQDATE)):
        if isExist == 0:
            status = 1
            break;
        else: 
            #附近有地震, 附近的地震都绘制过吗? 都没有绘制过, 就要画status = 1 一旦有绘制过的, 那就不画
            isPaint = isPainted4NearEq(EQDATE, selEQIndex, ind, d)
            if isPaint == 1: #都没有绘制过
               status = 1
               break;
    return status

#是否有近处的地震
def ExistNearEq(EQDATE, selEQIndex, d):
    status = 0 # 表示附近没有地震
    iEQDATE= EQDATE[selEQIndex]
    for i in range(len(EQDATE)):
        if i != selEQIndex and np.abs(iEQDATE - EQDATE[i]) < d:
           status = 1
           break;
    return status
#判断是否所有近处的地震都没有绘制过 status = 1
#判断是否所有近处的地震都有绘制过的 status = 0
def isPainted4NearEq(EQDATE, selEQIndex, ind, d):
    status = 1 #都没有绘制过;
    iEQDATE= EQDATE[selEQIndex]
    for i in range(len(EQDATE)):
        if i != selEQIndex and np.abs(iEQDATE - EQDATE[i]) < d and ind[i] == 1:
           status = 0
           break
    return status



def paintcomparefig4raw2pre(ouRoot, ouFile, RDDATE, RDDATA, PDDATE, PDDATA, yuzhi4up = None, yuzhi4dn = None):
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(RDDATE.T)
    fig, axs = plt.subplots(2, 1, figsize=(10, 4))
    axs[0].plot(RDDATE.T, RDDATA.T, marker = 'o',  linestyle = '-', color='blue', label="原始值", markersize = 0.5, linewidth = 0.2)
    axs[0].set_xlim(sdate, edate)
    axs[0].set_xticks(xtick4Index, xtick4Label)
    axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    axs[0].tick_params(axis='x', labelsize = 8)
    axs[0].tick_params(axis='y', labelsize = 8)
    axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[0].xaxis.set_minor_formatter(NullFormatter())

    axs[1].plot(PDDATE.T, PDDATA.T, marker = 'o',  linestyle = '-', color='blue', label="处理值", markersize = 0.5, linewidth = 0.2)
    if yuzhi4up != None:
       axs[1].plot([sdate, edate], [yuzhi4up, yuzhi4up], linestyle = '-', color='red', label="处理值", linewidth = 0.2)
    if yuzhi4dn != None:
       axs[1].plot([sdate, edate], [yuzhi4dn, yuzhi4dn], linestyle = '-', color='red', label="处理值", linewidth = 0.2)
    axs[1].set_xlim(sdate, edate)
    axs[1].set_xticks(xtick4Index, xtick4Label)
    #axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 3)
    axs[1].tick_params(axis='x', labelsize = 8)
    axs[1].tick_params(axis='y', labelsize = 8)
    axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[1].xaxis.set_minor_formatter(NullFormatter())
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print('输出文件:' + ouRoot + ouFile)
#绘制两条曲线，第一条为速率值；
#第二条为日均值;
def paintdisp4predictlongday4speed(DDATE, DDATA, DERRO, MDATE, MVALUE, TimeSpan, 
                             ouRoot, ouFile, 
                             yuzhi, 
                             NO, EQDATE, EQRANK, EQDIS, EQPOS, 
                             selrank, days4span, stnm):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''

    indices = MVALUE > yuzhi

    AB2V    = MVALUE[indices]
    AB2D    = MDATE[indices]
    #地震目录
    #print('地震目录')
    #print(EQDATE)
    al_bz4d, al_bz4v = py_tool.outputabnormaltimeserial(EQDATE, AB2D, AB2V, days4span)
    #print('异常线段')
    #print(len(al_bz4d))
    eq_bz4d, eq_bz4v = py_tool.outputcorrectedeq(EQDATE, AB2D, EQDIS, days4span)
    pr_bz4d, pr_bz4v = py_tool.outputpredictline(EQDATE, AB2D, AB2V, days4span)

    #还需要画一条近期预测曲线

    #还需要画一条近期预测曲线
    
    #print(TimeSpan[0])
    label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    #print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(DDATE.T)
    #up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    #up4D, lo4D                             = calYLim4Auto(DDATA, 0.0001)
    
    #print(xtick4Index)
    #print(xtick4Label)
    #fig, axs = plt.subplots(1, 1, figsize=(10, 4))
    #fig = plt.figure(figsize=(8, 2))
    #axs = fig.add_subplot(111)
    
    fig, axs = plt.subplots(2, 1, figsize=(10, 3))
    #axs[0].plot(TDATE, TDATA, color='blue',label="窗值1024s", linewidth = 0.2)
    #axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    #axs[0].set_xticks(xtick4Index, xtick4Label)
    #axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    #axs[0].tick_params(axis='x', labelsize = 8)
    #axs[0].tick_params(axis='y', labelsize = 8)
    #axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    #axs[0].xaxis.set_minor_formatter(NullFormatter())
    #axs.plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[0].plot(MDATE.T, MVALUE.T, marker = 'o', linestyle = '-', color='rosybrown', label=label4winL, markersize = 1, linewidth = 0.5)
    axs[0].plot([sdate, edate], [yuzhi, yuzhi], linestyle = '--', color='blue', linewidth = 0.5, label="阈值")
    axs[0].plot(AB2D, AB2V, marker = 'o', markersize = 2.5, linestyle = 'none', label="异常", markerfacecolor = 'none', markeredgecolor = 'red')
    axs[0].plot(al_bz4d, al_bz4v, linestyle = '-', color='red', label="持续时间", linewidth = 1)
    axs[0].plot(pr_bz4d, pr_bz4v, linestyle = '--', color='yellow', linewidth = 1, label = "预测")
    axs[0].set_ylabel(f"{stnm}" + '视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[0].set_xlabel('日期/yyyy', fontsize=8)
    axs[0].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs[0].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[0].xaxis.set_minor_formatter(NullFormatter())
    #axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.1, -0.29), ncol = 6)
    axs[0].legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.005, 1), ncol = 6)
    axs[0].tick_params(axis='x', labelsize = 8)
    axs[0].tick_params(axis='y', labelsize = 8)

    axs2 = axs[0].twinx()
    axs2.plot(EQDATE,  EQDIS,   color='black', label = '地震>{:.1f}'.format(selrank), marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'yellow')
    axs2.plot(eq_bz4d, eq_bz4v, color='black', label = '报准', marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'red')
    axs2.set_ylabel('震中距/km', fontsize= 8)
    axs2.tick_params(axis='y', labelcolor = 'red')
    disMin = 0
    disMax = np.max(EQDIS)
    disMax = int(disMax/100)
    disMax = (disMax + 2) * 100
    seq      = np.linspace(disMin, disMax, num = 5)
    axs2.set_yticks(seq)
    axs2.tick_params(axis='y', labelsize = 8)
    axs2.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    axs2.legend(prop={'size':6}, loc= 'lower right', bbox_to_anchor = (1.009, 1), ncol = 2)
    

    #判断这个时间序列能绘制多少个地震;
    EQRANK, EQDATE, EQPOS, EQDIS, ind = calpaintEqLabelIndex(sdate, edate, EQDATE, EQRANK, EQPOS, EQDIS)
    #相邻地震不足1个月的只绘制其中一个
    for i in range(len(EQDATE)):
        eq_date = EQDATE
        eq_valu = EQDIS
        if ind[i] == 1:
            axs2.annotate(
                "{}{:.1f}".format(EQPOS[i], EQRANK[i]),
                xy = (eq_date[i], eq_valu[i]),
                xytext = (3, 1),
                textcoords = 'offset points',
                arrowprops = None,
                fontsize   = 6,
                color='blue',
                bbox = None
            )
    axs[1].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[1].set_ylabel('视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[1].set_xlabel('日期/yyyy', fontsize=6)
    axs[1].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs[1].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[1].xaxis.set_minor_formatter(NullFormatter())
    axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    axs[1].tick_params(axis='x', labelsize = 8)
    axs[1].tick_params(axis='y', labelsize = 8)
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)
    #绘图后应该输出异常时间;
    #对预报的检验, 异常次数, 报准次数, 虚报次数
    #异常结构: abnormal_date, ab_value
    #报准结构: abnormal_corrected_date, abnormal_corrected_datd
    #虚报结构: abnormal_wrong_date, abnormal_wrong_data
    #对预报的地震, 地震次数, 报准次数, 漏报次数
    #总的地震结构: earthquake_date, earthquake_rank, earthquake_pos, earthquake_dis
    #报准地震结构: earthquake_corrected_date, earthquake_corrected_rank, earthquake_corrected_pos, earthquake_corrected_dis
    #报准地震结构: earthquake_wrong_date, earthquake_wrong_rank, earthquake_wrong_pos, earthquake_wrong_dis
    #R值 =  len(earthquake_corrected_date) / len(earthquake_date) - TimeSpan * len(abnormal_date) / (DDATE[-1] - DDATE[0])
    #R0值=  一个函数
    #需要写一段话;
def paintdisp4predictlongday4year(DDATE,   DDATA,     MDATE,  MVALUE, 
                                  ouRoot,  ouFile, 
                                  yuzhi, 
                                  NO,      EQDATE,    EQRANK, EQDIS, EQPOS, 
                                  selrank, days4span, stnm):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''
    
    indices = MVALUE[0, :] > yuzhi
    #print(indices)
    #print(MDATE)
    #print(MVALUE)
    AB2V    = MVALUE[0, indices]
    AB2D    = MDATE[0, indices]
    #地震目录
    #print('地震目录')
    #print(EQDATE)
    al_bz4d, al_bz4v = py_tool.outputabnormaltimeserial(EQDATE, AB2D, AB2V, days4span)
    #print('异常线段')
    #print(len(al_bz4d))
    eq_bz4d, eq_bz4v = py_tool.outputcorrectedeq(EQDATE, AB2D, EQDIS, days4span)
    pr_bz4d, pr_bz4v = py_tool.outputpredictline(EQDATE, AB2D, AB2V, days4span)

    #还需要画一条近期预测曲线

    #还需要画一条近期预测曲线
    
    #print(TimeSpan[0])
    #label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    #print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(DDATE.T)
    #up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    up4M, lo4M                             = calYLim4Auto2MaxMin(MVALUE[0])
    
    #print(xtick4Index)
    #print(xtick4Label)
    #fig, axs = plt.subplots(1, 1, figsize=(10, 4))
    fig = plt.figure(figsize=(8, 2))
    axs = fig.add_subplot(111)
    #fig, axs = plt.subplots(2, 1, figsize=(10, 4))
    #axs[0].plot(TDATE, TDATA, color='blue',label="窗值1024s", linewidth = 0.2)
    #axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    #axs[0].set_xticks(xtick4Index, xtick4Label)
    #axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    #axs[0].tick_params(axis='x', labelsize = 8)
    #axs[0].tick_params(axis='y', labelsize = 8)
    #axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    #axs[0].xaxis.set_minor_formatter(NullFormatter())
    axs.plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2, alpha = 0.3)
    axs.plot(MDATE.T, MVALUE[0].T, linestyle = '-', color='blue', label="模态1", linewidth = 0.8)
    axs.plot([sdate, edate], [yuzhi, yuzhi], linestyle = '--', color='red', linewidth = 0.5, label="阈值")
    
    x = MDATE[0].T
    y1= MVALUE[0].T
    y2= np.full(len(y1), yuzhi)
    
    axs.fill_between(x, y1, y2, where = (y1 >= y2), color = 'red', alpha = 0.2, label = '异常')
    axs.plot(AB2D, AB2V, marker = 'o', markersize = 0.2, linestyle = 'none', label="异常", markerfacecolor = 'none', markeredgecolor = 'red')
    axs.plot(al_bz4d, al_bz4v, linestyle = '-', color='red', label="持续时间", linewidth = 0.4, alpha = 0.3)
    axs.plot(pr_bz4d, pr_bz4v, linestyle = '--', color='yellow', linewidth = 1, label = "预测")
    axs.set_ylabel(f"{stnm}" + '视垂直位移/$10^{-6}m$', fontsize= 8)
    axs.set_xlabel('日期/yyyy', fontsize=8)
    axs.set_xlim(sdate, edate)
    axs.set_ylim(lo4M, up4M)
    axs.set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs.xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs.xaxis.set_minor_formatter(NullFormatter())
    #axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.1, -0.29), ncol = 6)
    axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.005, 1), ncol = 7)
    axs.tick_params(axis='x', labelsize = 8)
    axs.tick_params(axis='y', labelsize = 8)

    axs2 = axs.twinx()
    axs2.plot(EQDATE,  EQDIS,   color='black', label = '地震>{:.1f}'.format(selrank), marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'yellow')
    axs2.plot(eq_bz4d, eq_bz4v, color='black', label = '报准', marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'red')
    axs2.set_ylabel('震中距/km', fontsize= 8)
    axs2.tick_params(axis='y', labelcolor = 'red')
    disMin = 0
    disMax = np.max(EQDIS)
    disMax = int(disMax/100)
    disMax = (disMax + 2) * 100
    seq      = np.linspace(disMin, disMax, num = 5)
    axs2.set_yticks(seq)
    axs2.tick_params(axis='y', labelsize = 8)
    axs2.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    axs2.legend(prop={'size':6}, loc= 'lower right', bbox_to_anchor = (1.009, 1), ncol = 2)
    

    #判断这个时间序列能绘制多少个地震;
    EQRANK, EQDATE, EQPOS, EQDIS, ind = calpaintEqLabelIndex(sdate, edate, EQDATE, EQRANK, EQPOS, EQDIS)
    #相邻地震不足1个月的只绘制其中一个
    for i in range(len(EQDATE)):
        eq_date = EQDATE
        eq_valu = EQDIS
        if ind[i] == 1:
            axs2.annotate(
                "{}{:.1f}".format(EQPOS[i], EQRANK[i]),
                xy = (eq_date[i], eq_valu[i]),
                xytext = (3, 1),
                textcoords = 'offset points',
                arrowprops = None,
                fontsize   = 6,
                color='blue',
                bbox = None
            )
    #axs[1].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2)
    #axs[1].set_ylabel('视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[1].set_xlabel('日期/yyyy', fontsize=6)
    #axs[1].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    #axs[1].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    #axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    #axs[1].xaxis.set_minor_formatter(NullFormatter())
    #axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    #axs[1].tick_params(axis='x', labelsize = 8)
    #axs[1].tick_params(axis='y', labelsize = 8)
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)

def paintdisp4predictlongday4year4lo(DDATE,   DDATA,     MDATE,  MVALUE, 
                                  ouRoot,  ouFile, 
                                  yuzhi, 
                                  NO,      EQDATE,    EQRANK, EQDIS, EQPOS, 
                                  selrank, days4span, stnm):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''
    
    indices = MVALUE[0, :] < yuzhi
    #print(indices)
    #print(MDATE)
    #print(MVALUE)
    AB2V    = MVALUE[0, indices]
    AB2D    = MDATE[0, indices]
    #地震目录
    #print('地震目录')
    #print(EQDATE)
    al_bz4d, al_bz4v = py_tool.outputabnormaltimeserial(EQDATE, AB2D, AB2V, days4span)
    #print('测试异常线段')
    #print(al_bz4d)
    #print(al_bz4v)
    #print('测试异常线段')
    #print('异常线段')
    #print(len(al_bz4d))
    eq_bz4d, eq_bz4v = py_tool.outputcorrectedeq(EQDATE, AB2D, EQDIS, days4span)
    pr_bz4d, pr_bz4v = py_tool.outputpredictline(EQDATE, AB2D, AB2V, days4span)

    #还需要画一条近期预测曲线

    #还需要画一条近期预测曲线
    
    #print(TimeSpan[0])
    #label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    #print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(DDATE.T)
    #up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    up4M, lo4M                             = calYLim4Auto2MaxMin(MVALUE[0])
    
    #print(xtick4Index)
    #print(xtick4Label)
    #fig, axs = plt.subplots(1, 1, figsize=(10, 4))
    fig = plt.figure(figsize=(8, 2))
    axs = fig.add_subplot(111)
    #fig, axs = plt.subplots(2, 1, figsize=(10, 4))
    #axs[0].plot(TDATE, TDATA, color='blue',label="窗值1024s", linewidth = 0.2)
    #axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    #axs[0].set_xticks(xtick4Index, xtick4Label)
    #axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    #axs[0].tick_params(axis='x', labelsize = 8)
    #axs[0].tick_params(axis='y', labelsize = 8)
    #axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    #axs[0].xaxis.set_minor_formatter(NullFormatter())
    axs.plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2, alpha = 0.3)
    axs.plot(MDATE.T, MVALUE[0].T, linestyle = '-', color='blue', label="模态1", linewidth = 0.8)
    axs.plot([sdate, edate], [yuzhi, yuzhi], linestyle = '--', color='red', linewidth = 0.5, label="阈值")
    
    x = MDATE[0].T
    y1= MVALUE[0].T
    y2= np.full(len(y1), yuzhi)
    
    axs.fill_between(x, y1, y2, where = (y1 <= y2), color = 'red', alpha = 0.2, label = '异常')
    axs.plot(AB2D, AB2V, marker = 'o', markersize = 0.2, linestyle = 'none', label="异常", markerfacecolor = 'none', markeredgecolor = 'red')
    axs.plot(al_bz4d, al_bz4v, linestyle = '-',  color='red',    label="持续时间", linewidth = 0.4, alpha = 0.3)
    axs.plot(pr_bz4d, pr_bz4v, linestyle = '--', color='yellow', linewidth = 1,   label = "预测")
    axs.set_ylabel(f"{stnm}" + '视垂直位移/$10^{-6}m$', fontsize= 8)
    axs.set_xlabel('日期/yyyy', fontsize=8)
    axs.set_xlim(sdate, edate)
    axs.set_ylim(lo4M, up4M)
    axs.set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs.xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs.xaxis.set_minor_formatter(NullFormatter())
    #axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.1, -0.29), ncol = 6)
    axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.005, 1), ncol = 7)
    axs.tick_params(axis='x', labelsize = 8)
    axs.tick_params(axis='y', labelsize = 8)

    axs2 = axs.twinx()
    axs2.plot(EQDATE,  EQDIS,   color='black', label = '地震>{:.1f}'.format(selrank), marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'yellow')
    axs2.plot(eq_bz4d, eq_bz4v, color='black', label = '报准', marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'red')
    axs2.set_ylabel('震中距/km', fontsize= 8)
    axs2.tick_params(axis='y', labelcolor = 'red')
    disMin = 0
    disMax = np.max(EQDIS)
    disMax = int(disMax/100)
    disMax = (disMax + 2) * 100
    seq      = np.linspace(disMin, disMax, num = 5)
    axs2.set_yticks(seq)
    axs2.tick_params(axis='y', labelsize = 8)
    axs2.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    axs2.legend(prop={'size':6}, loc= 'lower right', bbox_to_anchor = (1.009, 1), ncol = 2)
    

    #判断这个时间序列能绘制多少个地震;
    EQRANK, EQDATE, EQPOS, EQDIS, ind = calpaintEqLabelIndex(sdate, edate, EQDATE, EQRANK, EQPOS, EQDIS)
    #相邻地震不足1个月的只绘制其中一个
    for i in range(len(EQDATE)):
        eq_date = EQDATE
        eq_valu = EQDIS
        if ind[i] == 1:
            axs2.annotate(
                "{}{:.1f}".format(EQPOS[i], EQRANK[i]),
                xy = (eq_date[i], eq_valu[i]),
                xytext = (3, 1),
                textcoords = 'offset points',
                arrowprops = None,
                fontsize   = 6,
                color='blue',
                bbox = None
            )
    #axs[1].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2)
    #axs[1].set_ylabel('视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[1].set_xlabel('日期/yyyy', fontsize=6)
    #axs[1].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    #axs[1].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    #axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    #axs[1].xaxis.set_minor_formatter(NullFormatter())
    #axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    #axs[1].tick_params(axis='x', labelsize = 8)
    #axs[1].tick_params(axis='y', labelsize = 8)
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)

def paintdisp4predictlongday(DDATE, DDATA, DERRO, MDATE, MVALUE, TimeSpan, 
                             ouRoot, ouFile, 
                             yuzhi, 
                             NO, EQDATE, EQRANK, EQDIS, EQPOS, 
                             selrank, days4span, stnm):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''

    indices = MVALUE > yuzhi

    AB2V    = MVALUE[indices]
    AB2D    = MDATE[indices]
    #地震目录
    #print('地震目录')
    #print(EQDATE)
    al_bz4d, al_bz4v = py_tool.outputabnormaltimeserial(EQDATE, AB2D, AB2V, days4span)
    #print('异常线段')
    #print(len(al_bz4d))
    eq_bz4d, eq_bz4v = py_tool.outputcorrectedeq(EQDATE, AB2D, EQDIS, days4span)
    pr_bz4d, pr_bz4v = py_tool.outputpredictline(EQDATE, AB2D, AB2V, days4span)

    #还需要画一条近期预测曲线

    #还需要画一条近期预测曲线
    
    #print(TimeSpan[0])
    label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    #print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(DDATE.T)
    #up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    #up4D, lo4D                             = calYLim4Auto(DDATA, 0.0001)
    
    #print(xtick4Index)
    #print(xtick4Label)
    #fig, axs = plt.subplots(1, 1, figsize=(10, 4))
    fig = plt.figure(figsize=(8, 2))
    axs = fig.add_subplot(111)
    #axs[0].plot(TDATE, TDATA, color='blue',label="窗值1024s", linewidth = 0.2)
    #axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    #axs[0].set_xticks(xtick4Index, xtick4Label)
    #axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    #axs[0].tick_params(axis='x', labelsize = 8)
    #axs[0].tick_params(axis='y', labelsize = 8)
    #axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    #axs[0].xaxis.set_minor_formatter(NullFormatter())
    axs.plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs.plot(MDATE.T, MVALUE.T, marker = 'o', linestyle = '-', color='rosybrown', label=label4winL, markersize = 1, linewidth = 0.5)
    axs.plot([sdate, edate], [yuzhi, yuzhi], linestyle = '--', color='blue', linewidth = 0.5, label="阈值")
    axs.plot(AB2D, AB2V, marker = 'o', markersize = 2.5, linestyle = 'none', label="异常", markerfacecolor = 'none', markeredgecolor = 'red')
    axs.plot(al_bz4d, al_bz4v, linestyle = '-', color='red', label="持续时间", linewidth = 1)
    axs.plot(pr_bz4d, pr_bz4v, linestyle = '--', color='yellow', linewidth = 1, label = "预测")
    axs.set_ylabel(f"{stnm}" + '视垂直位移/$10^{-6}m$', fontsize= 8)
    axs.set_xlabel('日期/yyyy', fontsize=8)
    axs.set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs.set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs.xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs.xaxis.set_minor_formatter(NullFormatter())
    #axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.1, -0.29), ncol = 6)
    axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.005, 1), ncol = 6)
    axs.tick_params(axis='x', labelsize = 8)
    axs.tick_params(axis='y', labelsize = 8)

    axs2 = axs.twinx()
    axs2.plot(EQDATE,  EQDIS,   color='black', label = '地震>{:.1f}'.format(selrank), marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'yellow')
    axs2.plot(eq_bz4d, eq_bz4v, color='black', label = '报准', marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'red')
    axs2.set_ylabel('震中距/km', fontsize= 8)
    axs2.tick_params(axis='y', labelcolor = 'red')
    disMin = 0
    disMax = np.max(EQDIS)
    disMax = int(disMax/100)
    disMax = (disMax + 1) * 100
    seq      = np.linspace(disMin, disMax, num = 5)
    axs2.set_yticks(seq)
    axs2.tick_params(axis='y', labelsize = 8)
    axs2.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    axs2.legend(prop={'size':6}, loc= 'lower right', bbox_to_anchor = (1.009, 1), ncol = 2)
    

    #判断这个时间序列能绘制多少个地震;
    EQRANK, EQDATE, EQPOS, EQDIS, ind = calpaintEqLabelIndex(sdate, edate, EQDATE, EQRANK, EQPOS, EQDIS)
    #相邻地震不足1个月的只绘制其中一个
    for i in range(len(EQDATE)):
        eq_date = EQDATE
        eq_valu = EQDIS
        if ind[i] == 1:
            axs2.annotate(
                "{}{:.1f}".format(EQPOS[i], EQRANK[i]),
                xy = (eq_date[i], eq_valu[i]),
                xytext = (3, 1),
                textcoords = 'offset points',
                arrowprops = None,
                fontsize   = 6,
                color='blue',
                bbox = None
            )
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)
    #绘图后应该输出异常时间;
    #对预报的检验, 异常次数, 报准次数, 虚报次数
    #异常结构: abnormal_date, ab_value
    #报准结构: abnormal_corrected_date, abnormal_corrected_datd
    #虚报结构: abnormal_wrong_date, abnormal_wrong_data
    #对预报的地震, 地震次数, 报准次数, 漏报次数
    #总的地震结构: earthquake_date, earthquake_rank, earthquake_pos, earthquake_dis
    #报准地震结构: earthquake_corrected_date, earthquake_corrected_rank, earthquake_corrected_pos, earthquake_corrected_dis
    #报准地震结构: earthquake_wrong_date, earthquake_wrong_rank, earthquake_wrong_pos, earthquake_wrong_dis
    #R值 =  len(earthquake_corrected_date) / len(earthquake_date) - TimeSpan * len(abnormal_date) / (DDATE[-1] - DDATE[0])
    #R0值=  一个函数
    #需要写一段话;

def paintdisp4splitlongday(TDATE, TDATA, DDATE, DDATA, DERRO, MDATE, MVALUE, TimeSpan, ouRoot, ouFile):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''
    
    #print(TimeSpan[0])
    label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(TDATE)
    up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    up4D, lo4D                             = calYLim4Auto(DDATA, 0.0001)
    
    print(xtick4Index)
    print(xtick4Label)
    fig, axs = plt.subplots(2, 1, figsize=(10, 4))

    axs[0].plot(TDATE, TDATA, color='blue',label="窗值1024s", linewidth = 0.2)
    axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    axs[0].set_xticks(xtick4Index, xtick4Label)
    axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    axs[0].tick_params(axis='x', labelsize = 8)
    axs[0].tick_params(axis='y', labelsize = 8)
    axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[0].xaxis.set_minor_formatter(NullFormatter())
    axs[1].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='gray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[1].plot(MDATE.T, MVALUE.T, marker = 'o', linestyle = '-', color='red', label=label4winL, markersize = 0.5, linewidth = 0.2)
    axs[1].set_ylabel('视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[1].set_xlabel('日期/yyyy', fontsize=6)
    axs[1].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs[1].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[1].xaxis.set_minor_formatter(NullFormatter())
    axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    axs[1].tick_params(axis='x', labelsize = 8)
    axs[1].tick_params(axis='y', labelsize = 8)
    
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)

def calpaintEqLabelIndex(SDATE, EDATE, EQDATE, EQRANK, EQPOS, EQDIS):
    N      = EDATE - SDATE
    Y      = int(N/365)
    Y3     = Y * 3 #有多少个地震都只绘制这么多地震;
    M      = len(EQRANK)#震都只绘制这么多地震;
    M      = len(EQRANK)
    d      = 200
    ind    = np.full(M, 0)    #如果序列为1,表示要绘制, 否者不绘制
    sorted_indices = np.argsort(-EQRANK)
    EQRANK         = EQRANK[sorted_indices]
    EQDATE         = EQDATE[sorted_indices]
    EQPOS          = EQPOS[sorted_indices]
    EQDIS          = EQDIS[sorted_indices]

    for i in range(len(ind)):
        if isANNO(EQDATE, i, d, ind) == 1 and i <= Y3:
            ind[i] = 1
    return EQRANK, EQDATE, EQPOS, EQDIS, ind

#判断列表中是否有比该地震更大的地震
def isANNO(EQDATE, selEQIndex, d, ind):
    status = 0 # 不画
    iEQDATE = EQDATE[selEQIndex]
    isExist = ExistNearEq(EQDATE, selEQIndex, d); # isExist = 0 附近没有地震; isExist = 1 附近有地震;
    for i in range(len(EQDATE)):
        if isExist == 0:
            status = 1
            break;
        else: 
            #附近有地震, 附近的地震都绘制过吗? 都没有绘制过, 就要画status = 1 一旦有绘制过的, 那就不画
            isPaint = isPainted4NearEq(EQDATE, selEQIndex, ind, d)
            if isPaint == 1: #都没有绘制过
               status = 1
               break;
    return status

#是否有近处的地震
def ExistNearEq(EQDATE, selEQIndex, d):
    status = 0 # 表示附近没有地震
    iEQDATE= EQDATE[selEQIndex]
    for i in range(len(EQDATE)):
        if i != selEQIndex and np.abs(iEQDATE - EQDATE[i]) < d:
           status = 1
           break;
    return status
#判断是否所有近处的地震都没有绘制过 status = 1
#判断是否所有近处的地震都有绘制过的 status = 0
def isPainted4NearEq(EQDATE, selEQIndex, ind, d):
    status = 1 #都没有绘制过;
    iEQDATE= EQDATE[selEQIndex]
    for i in range(len(EQDATE)):
        if i != selEQIndex and np.abs(iEQDATE - EQDATE[i]) < d and ind[i] == 1:
           status = 0
           break
    return status



def paintcomparefig4raw2pre(ouRoot, ouFile, RDDATE, RDDATA, PDDATE, PDDATA, yuzhi4up = None, yuzhi4dn = None):
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(RDDATE.T)
    fig, axs = plt.subplots(2, 1, figsize=(10, 4))
    axs[0].plot(RDDATE.T, RDDATA.T, marker = 'o',  linestyle = '-', color='blue', label="原始值", markersize = 0.5, linewidth = 0.2)
    axs[0].set_xlim(sdate, edate)
    axs[0].set_xticks(xtick4Index, xtick4Label)
    axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    axs[0].tick_params(axis='x', labelsize = 8)
    axs[0].tick_params(axis='y', labelsize = 8)
    axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[0].xaxis.set_minor_formatter(NullFormatter())

    axs[1].plot(PDDATE.T, PDDATA.T, marker = 'o',  linestyle = '-', color='blue', label="处理值", markersize = 0.5, linewidth = 0.2)
    if yuzhi4up != None:
       axs[1].plot([sdate, edate], [yuzhi4up, yuzhi4up], linestyle = '-', color='red', label="处理值", linewidth = 0.2)
    if yuzhi4dn != None:
       axs[1].plot([sdate, edate], [yuzhi4dn, yuzhi4dn], linestyle = '-', color='red', label="处理值", linewidth = 0.2)
    axs[1].set_xlim(sdate, edate)
    axs[1].set_xticks(xtick4Index, xtick4Label)
    #axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 3)
    axs[1].tick_params(axis='x', labelsize = 8)
    axs[1].tick_params(axis='y', labelsize = 8)
    axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[1].xaxis.set_minor_formatter(NullFormatter())
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print('输出文件:' + ouRoot + ouFile)
#绘制两条曲线，第一条为速率值；
#第二条为日均值;
def paintdisp4predictlongday4speed(DDATE, DDATA, DERRO, MDATE, MVALUE, TimeSpan, 
                             ouRoot, ouFile, 
                             yuzhi, 
                             NO, EQDATE, EQRANK, EQDIS, EQPOS, 
                             selrank, days4span, stnm):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''

    indices = MVALUE > yuzhi

    AB2V    = MVALUE[indices]
    AB2D    = MDATE[indices]
    #地震目录
    #print('地震目录')
    #print(EQDATE)
    al_bz4d, al_bz4v = py_tool.outputabnormaltimeserial(EQDATE, AB2D, AB2V, days4span)
    #print('异常线段')
    #print(len(al_bz4d))
    eq_bz4d, eq_bz4v = py_tool.outputcorrectedeq(EQDATE, AB2D, EQDIS, days4span)
    pr_bz4d, pr_bz4v = py_tool.outputpredictline(EQDATE, AB2D, AB2V, days4span)

    #还需要画一条近期预测曲线

    #还需要画一条近期预测曲线
    
    #print(TimeSpan[0])
    label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    #print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(DDATE.T)
    #up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    #up4D, lo4D                             = calYLim4Auto(DDATA, 0.0001)
    
    #print(xtick4Index)
    #print(xtick4Label)
    #fig, axs = plt.subplots(1, 1, figsize=(10, 4))
    #fig = plt.figure(figsize=(8, 2))
    #axs = fig.add_subplot(111)
    fig, axs = plt.subplots(2, 1, figsize=(10, 4))
    #axs[0].plot(TDATE, TDATA, color='blue',label="窗值1024s", linewidth = 0.2)
    #axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    #axs[0].set_xticks(xtick4Index, xtick4Label)
    #axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    #axs[0].tick_params(axis='x', labelsize = 8)
    #axs[0].tick_params(axis='y', labelsize = 8)
    #axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    #axs[0].xaxis.set_minor_formatter(NullFormatter())
    #axs.plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[0].plot(MDATE.T, MVALUE.T, marker = 'o', linestyle = '-', color='rosybrown', label=label4winL, markersize = 1, linewidth = 0.5)
    axs[0].plot([sdate, edate], [yuzhi, yuzhi], linestyle = '--', color='blue', linewidth = 0.5, label="阈值")
    axs[0].plot(AB2D, AB2V, marker = 'o', markersize = 2.5, linestyle = 'none', label="异常", markerfacecolor = 'none', markeredgecolor = 'red')
    axs[0].plot(al_bz4d, al_bz4v, linestyle = '-', color='red', label="持续时间", linewidth = 1)
    axs[0].plot(pr_bz4d, pr_bz4v, linestyle = '--', color='yellow', linewidth = 1, label = "预测")
    axs[0].set_ylabel(f"{stnm}" + '视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[0].set_xlabel('日期/yyyy', fontsize=8)
    axs[0].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs[0].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[0].xaxis.set_minor_formatter(NullFormatter())
    #axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.1, -0.29), ncol = 6)
    axs[0].legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.005, 1), ncol = 6)
    axs[0].tick_params(axis='x', labelsize = 8)
    axs[0].tick_params(axis='y', labelsize = 8)

    axs2 = axs[0].twinx()
    axs2.plot(EQDATE,  EQDIS,   color='black', label = '地震>{:.1f}'.format(selrank), marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'yellow')
    axs2.plot(eq_bz4d, eq_bz4v, color='black', label = '报准', marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'red')
    axs2.set_ylabel('震中距/km', fontsize= 8)
    axs2.tick_params(axis='y', labelcolor = 'red')
    disMin = 0
    disMax = np.max(EQDIS)
    disMax = int(disMax/100)
    disMax = (disMax + 2) * 100
    seq      = np.linspace(disMin, disMax, num = 5)
    axs2.set_yticks(seq)
    axs2.tick_params(axis='y', labelsize = 8)
    axs2.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    axs2.legend(prop={'size':6}, loc= 'lower right', bbox_to_anchor = (1.009, 1), ncol = 2)
    

    #判断这个时间序列能绘制多少个地震;
    EQRANK, EQDATE, EQPOS, EQDIS, ind = calpaintEqLabelIndex(sdate, edate, EQDATE, EQRANK, EQPOS, EQDIS)
    #相邻地震不足1个月的只绘制其中一个
    for i in range(len(EQDATE)):
        eq_date = EQDATE
        eq_valu = EQDIS
        if ind[i] == 1:
            axs2.annotate(
                "{}{:.1f}".format(EQPOS[i], EQRANK[i]),
                xy = (eq_date[i], eq_valu[i]),
                xytext = (3, 1),
                textcoords = 'offset points',
                arrowprops = None,
                fontsize   = 6,
                color='blue',
                bbox = None
            )
    axs[1].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[1].set_ylabel('视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[1].set_xlabel('日期/yyyy', fontsize=6)
    axs[1].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs[1].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[1].xaxis.set_minor_formatter(NullFormatter())
    axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    axs[1].tick_params(axis='x', labelsize = 8)
    axs[1].tick_params(axis='y', labelsize = 8)
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)
    #绘图后应该输出异常时间;
    #对预报的检验, 异常次数, 报准次数, 虚报次数
    #异常结构: abnormal_date, ab_value
    #报准结构: abnormal_corrected_date, abnormal_corrected_datd
    #虚报结构: abnormal_wrong_date, abnormal_wrong_data
    #对预报的地震, 地震次数, 报准次数, 漏报次数
    #总的地震结构: earthquake_date, earthquake_rank, earthquake_pos, earthquake_dis
    #报准地震结构: earthquake_corrected_date, earthquake_corrected_rank, earthquake_corrected_pos, earthquake_corrected_dis
    #报准地震结构: earthquake_wrong_date, earthquake_wrong_rank, earthquake_wrong_pos, earthquake_wrong_dis
    #R值 =  len(earthquake_corrected_date) / len(earthquake_date) - TimeSpan * len(abnormal_date) / (DDATE[-1] - DDATE[0])
    #R0值=  一个函数
    #需要写一段话;
def paintdisp4predictlongday(DDATE, DDATA, DERRO, MDATE, MVALUE, TimeSpan, 
                             ouRoot, ouFile, 
                             yuzhi, 
                             NO, EQDATE, EQRANK, EQDIS, EQPOS, 
                             selrank, days4span, stnm):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''

    indices = MVALUE > yuzhi

    AB2V    = MVALUE[indices]
    AB2D    = MDATE[indices]
    #地震目录
    #print('地震目录')
    #print(EQDATE)
    al_bz4d, al_bz4v = py_tool.outputabnormaltimeserial(EQDATE, AB2D, AB2V, days4span)
    #print('异常线段')
    #print(len(al_bz4d))
    eq_bz4d, eq_bz4v = py_tool.outputcorrectedeq(EQDATE, AB2D, EQDIS, days4span)
    pr_bz4d, pr_bz4v = py_tool.outputpredictline(EQDATE, AB2D, AB2V, days4span)

    #还需要画一条近期预测曲线

    #还需要画一条近期预测曲线
    
    #print(TimeSpan[0])
    label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    #print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(DDATE.T)
    #up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    #up4D, lo4D                             = calYLim4Auto(DDATA, 0.0001)
    
    #print(xtick4Index)
    #print(xtick4Label)
    #fig, axs = plt.subplots(1, 1, figsize=(10, 4))
    matplotlib.rc("font",family='SimSun')
    fig = plt.figure(figsize=(8, 2))
    axs = fig.add_subplot(111)
    #axs[0].plot(TDATE, TDATA, color='blue',label="窗值1024s", linewidth = 0.2)
    #axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    #axs[0].set_xticks(xtick4Index, xtick4Label)
    #axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    #axs[0].tick_params(axis='x', labelsize = 8)
    #axs[0].tick_params(axis='y', labelsize = 8)
    #axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    #axs[0].xaxis.set_minor_formatter(NullFormatter())
    axs.plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='lightgray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs.plot(MDATE.T, MVALUE.T, marker = 'o', linestyle = '-', color='rosybrown', label=label4winL, markersize = 1, linewidth = 0.5)
    axs.plot([sdate, edate], [yuzhi, yuzhi], linestyle = '--', color='blue', linewidth = 0.5, label="阈值")
    axs.plot(AB2D, AB2V, marker = 'o', markersize = 2.5, linestyle = 'none', label="异常", markerfacecolor = 'none', markeredgecolor = 'red')
    axs.plot(al_bz4d, al_bz4v, linestyle = '-', color='red', label="持续时间", linewidth = 1)
    axs.plot(pr_bz4d, pr_bz4v, linestyle = '--', color='yellow', linewidth = 1, label = "预测")
    axs.set_ylabel(f"{stnm}" + '视垂直位移/$10^{-6}m$', fontsize= 8)
    axs.set_xlabel('日期/yyyy', fontsize=8)
    axs.set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs.set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs.xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs.xaxis.set_minor_formatter(NullFormatter())
    #axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.1, -0.29), ncol = 6)
    axs.legend(prop={'size':6}, loc= 'lower left', bbox_to_anchor = (-0.005, 1), ncol = 6)
    axs.tick_params(axis='x', labelsize = 8)
    axs.tick_params(axis='y', labelsize = 8)

    axs2 = axs.twinx()
    axs2.plot(EQDATE,  EQDIS,   color='black', label = '地震>{:.1f}'.format(selrank), marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'yellow')
    axs2.plot(eq_bz4d, eq_bz4v, color='black', label = '报准', marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'red')
    axs2.set_ylabel('震中距/km', fontsize= 8)
    axs2.tick_params(axis='y', labelcolor = 'red')
    disMin = 0
    disMax = np.max(EQDIS)
    disMax = int(disMax/100)
    disMax = (disMax + 1) * 100
    seq      = np.linspace(disMin, disMax, num = 5)
    axs2.set_yticks(seq)
    axs2.tick_params(axis='y', labelsize = 8)
    axs2.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    axs2.legend(prop={'size':6}, loc= 'lower right', bbox_to_anchor = (1.009, 1), ncol = 2)
    

    #判断这个时间序列能绘制多少个地震;
    EQRANK, EQDATE, EQPOS, EQDIS, ind = calpaintEqLabelIndex(sdate, edate, EQDATE, EQRANK, EQPOS, EQDIS)
    #相邻地震不足1个月的只绘制其中一个
    for i in range(len(EQDATE)):
        eq_date = EQDATE
        eq_valu = EQDIS
        if ind[i] == 1:
            axs2.annotate(
                "{}{:.1f}".format(EQPOS[i], EQRANK[i]),
                xy = (eq_date[i], eq_valu[i]),
                xytext = (3, 1),
                textcoords = 'offset points',
                arrowprops = None,
                fontsize   = 6,
                color='blue',
                bbox = None
            )
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)
    #绘图后应该输出异常时间;
    #对预报的检验, 异常次数, 报准次数, 虚报次数
    #异常结构: abnormal_date, ab_value
    #报准结构: abnormal_corrected_date, abnormal_corrected_datd
    #虚报结构: abnormal_wrong_date, abnormal_wrong_data
    #对预报的地震, 地震次数, 报准次数, 漏报次数
    #总的地震结构: earthquake_date, earthquake_rank, earthquake_pos, earthquake_dis
    #报准地震结构: earthquake_corrected_date, earthquake_corrected_rank, earthquake_corrected_pos, earthquake_corrected_dis
    #报准地震结构: earthquake_wrong_date, earthquake_wrong_rank, earthquake_wrong_pos, earthquake_wrong_dis
    #R值 =  len(earthquake_corrected_date) / len(earthquake_date) - TimeSpan * len(abnormal_date) / (DDATE[-1] - DDATE[0])
    #R0值=  一个函数
    #需要写一段话;

def paintdisp4vmdlongday(TDATE, TDATA, DDATE, DDATA, ouRoot, ouFile):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''
    
    #print(TimeSpan[0])
    #label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    #print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(DDATE.T)
    up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    up4D, lo4D                             = calYLim4Auto(DDATA, 0.0001)
    
    print(xtick4Index)
    print(xtick4Label)
    fig, axs = plt.subplots(2, 1, figsize=(10, 4))
    row, col = TDATA.shape
    colorlist = ['red', 'blue', 'green', 'yellow', 'gray']
    for i in range(row):
        if i <= 3:
            #print(TDATE)
            #print(TDATA[i])
            axs[0].plot(TDATE[0], TDATA[i], color=colorlist[i],label="模态{:d}".format(i + 1 ), linewidth = 0.2)

    axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    axs[0].set_xticks(xtick4Index, xtick4Label)
    axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 5)
    axs[0].tick_params(axis='x', labelsize = 8)
    axs[0].tick_params(axis='y', labelsize = 8)
    axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[0].xaxis.set_minor_formatter(NullFormatter())
    axs[1].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='gray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[1].plot(TDATE.T, TDATA[0, :].T, linestyle = '-', color='red', label="模态0", linewidth = 0.8)
    axs[1].set_ylabel('视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[1].set_xlabel('日期/yyyy', fontsize=6)
    axs[1].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs[1].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[1].xaxis.set_minor_formatter(NullFormatter())
    axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    axs[1].tick_params(axis='x', labelsize = 8)
    axs[1].tick_params(axis='y', labelsize = 8)
    
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)

def paintdisp4splitlongday2speed(TDATE, TDATA, DDATE, DDATA, DERRO, MDATE, MVALUE, PDATE, PDATA, TimeSpan, ouRoot, ouFile):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''
    
    #print(TimeSpan[0])
    label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(TDATE)
    up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    up4D, lo4D                             = calYLim4Auto(DDATA, 0.0001)
    
    print(xtick4Index)
    print(xtick4Label)
    fig, axs = plt.subplots(3, 1, figsize=(10, 5))

    axs[0].plot(TDATE, TDATA, color='blue',label="窗值1024s", linewidth = 0.2)
    axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    axs[0].set_xticks(xtick4Index, xtick4Label)
    axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    axs[0].tick_params(axis='x', labelsize = 8)
    axs[0].tick_params(axis='y', labelsize = 8)
    axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[0].xaxis.set_minor_formatter(NullFormatter())
    axs[1].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='gray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[1].plot(MDATE.T, MVALUE.T, marker = 'o', linestyle = '-', color='red', label=label4winL, markersize = 0.5, linewidth = 0.2)
    axs[1].set_ylabel('视垂直位移/$10^{-6}m$', fontsize= 8)
    axs[1].set_xlabel('日期/yyyy', fontsize=6)
    axs[1].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs[1].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[1].xaxis.set_minor_formatter(NullFormatter())
    axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    axs[1].tick_params(axis='x', labelsize = 8)
    axs[1].tick_params(axis='y', labelsize = 8)


    #axs[1].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='gray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[2].plot(PDATE.T, PDATA.T, marker = 'o', linestyle = '-', color='red', label=label4winL, markersize = 0.5, linewidth = 0.2)
    axs[2].set_ylabel('速率/$10^{-8}m/M$', fontsize= 8)
    axs[2].set_xlabel('日期/yyyy', fontsize=6)
    axs[2].set_xlim(sdate, edate)
    axs[2].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[2].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[2].xaxis.set_minor_formatter(NullFormatter())
    axs[2].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    axs[2].tick_params(axis='x', labelsize = 8)
    axs[2].tick_params(axis='y', labelsize = 8)

    #双Y轴绘图

    #双Y轴绘图
    
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)

def paintdisp4splitlongday2speed4eq(TDATE, TDATA, 
    DDATE, DDATA, DERRO, 
    MDATE, MVALUE, 
    PDATE, PDATA, 
    TimeSpan, ouRoot, ouFile,
    NO,    EQDATE,   EQRANK,  EQDIS,  EQPOS, stnm, selrank):
    '''
    TDATE和DDATA是时间对应的浮点数;
    DS给出用pd的方法来绘制时间序列，我也用;
    '''
    
    #print(TimeSpan[0])
    label4winL = "窗长{}d".format(TimeSpan[0,0].tolist())
    print(label4winL)
    sdate, edate, xtick4Index, xtick4Label, xtick4Month = sdate2edate4timeserial4any(TDATE)
    up4T, lo4T                             = calYLim4Auto(TDATA, 0.0001)
    up4D, lo4D                             = calYLim4Auto(DDATA, 0.0001)
    
    print(xtick4Index)
    print(xtick4Label)
    fig, axs = plt.subplots(2, 1, figsize=(10, 3))

    #axs[0].plot(TDATE, TDATA, color='blue',label="窗值1024s", linewidth = 0.2)
    #axs[0].set_ylabel(r'视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[0].set_xlim(sdate, edate)
    #axs[0].set_ylim(lo4T, up4T)
    #axs[0].set_xticks(xtick4Index, xtick4Label)
    #axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 1)
    #axs[0].tick_params(axis='x', labelsize = 8)
    #axs[0].tick_params(axis='y', labelsize = 8)
    #axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    #axs[0].xaxis.set_minor_formatter(NullFormatter())
    axs[0].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='gray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[0].plot(MDATE.T, MVALUE.T, marker = 'o', linestyle = '-', color='red', label=label4winL, markersize = 0.5, linewidth = 0.2)
    axs[0].set_ylabel(stnm + '视垂直位移/$10^{-6}m$', fontsize= 8)
    #axs[0].set_xlabel('日期/yyyy', fontsize=6)
    axs[0].set_xlim(sdate, edate)
    #axs[1].set_ylim(lo4D, up4D)
    axs[0].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[0].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[0].xaxis.set_minor_formatter(NullFormatter())
    axs[0].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    axs[0].tick_params(axis='x', labelsize = 8)
    axs[0].tick_params(axis='y', labelsize = 8)
    axs2 = axs[0].twinx()
    axs2.plot(EQDATE,  EQDIS,   color='black', label = '地震>{:.1f}'.format(selrank), marker = 'o', linestyle = 'none', markersize = 6, markerfacecolor = 'red', markeredgecolor = 'red')
    #axs2.plot(eq_bz4d, eq_bz4v, color='black', label = '报准', marker = 'o', linestyle = 'none', markersize = 4, markerfacecolor = 'none', markeredgecolor = 'red')
    axs2.set_ylabel('震中距/km', fontsize= 8)
    axs2.tick_params(axis='y', labelcolor = 'red')
    disMin = 0
    disMax = np.max(EQDIS)
    disMax = int(disMax/100)
    disMax = (disMax + 2) * 100
    seq      = np.linspace(disMin, disMax, num = 5)
    axs2.set_yticks(seq)
    axs2.tick_params(axis='y', labelsize = 8)
    axs2.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    axs2.legend(prop={'size':6}, loc= 'upper right', ncol = 2)
    

    #判断这个时间序列能绘制多少个地震;
    EQRANK, EQDATE, EQPOS, EQDIS, ind = calpaintEqLabelIndex(sdate, edate, EQDATE, EQRANK, EQPOS, EQDIS)
    #相邻地震不足1个月的只绘制其中一个
    for i in range(len(EQDATE)):
        eq_date = EQDATE
        eq_valu = EQDIS
        if ind[i] == 1:
            axs2.annotate(
                "{}{:.1f}".format(EQPOS[i], EQRANK[i]),
                xy = (eq_date[i], eq_valu[i]),
                xytext = (3, 1),
                textcoords = 'offset points',
                arrowprops = None,
                fontsize   = 6,
                color='blue',
                bbox = None
            )




    #axs[1].plot(DDATE.T, DDATA.T, marker = 'o',  linestyle = '-', color='gray', label="日值", markersize = 0.5, linewidth = 0.2)
    axs[1].plot(PDATE.T, PDATA.T, marker = 'o', linestyle = '-', color='red', label=label4winL, markersize = 0.5, linewidth = 0.2)
    axs[1].set_ylabel('速率/$10^{-8}m/M$', fontsize= 8)
    axs[1].set_xlabel('日期/yyyy', fontsize=6)
    axs[1].set_xlim(sdate, edate)
    axs[1].set_xticks(xtick4Index, xtick4Label, fontsize= 8)
    axs[1].xaxis.set_minor_locator(FixedLocator(xtick4Month))
    axs[1].xaxis.set_minor_formatter(NullFormatter())
    axs[1].legend(prop={'size':6}, loc= 'upper left', ncol = 2)
    axs[1].tick_params(axis='x', labelsize = 8)
    axs[1].tick_params(axis='y', labelsize = 8)

    #双Y轴绘图

    #双Y轴绘图
    
    plt.savefig(ouRoot + ouFile, dpi=250, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)


def paintdisp4longday(TDATE, TDATA, DDATE, DDATA, DERRO, ouRoot, ouFile):
    plt.figure(figsize=(10,6))
    plt.plot(TDATE, TDATA, color='blue',label="位移")
    #plt.plot(DDATE, DDATA, marker = 'o', linestyle = '-', color = 'r', label = '日值', markersize = 1)
    plt.xlabel('Date[yy-mm-dd]', fontsize=10)
    plt.ylabel('Disp/um', fontsize=12)
    plt.title('Time Series Plot')
    sdate, edate, xtick4Index, xtick4Label = sdate2edate4timeserial(TDATE)
    plt.xticks(xtick4Index, xtick4Label)  # 设置X轴刻度位置和标签
    plt.ylim(0, 0.3)
    # 显示图表（可选）
    plt.savefig(ouRoot + ouFile, dpi=300, bbox_inches='tight')
    plt.show()
    print(ouRoot + ouFile)

def paintdisp4win2day(TDATE, TDATA, DDATE, DDATA, DERRO):
    ouRoot= '/home/mw/temp/'
    #matplotlib.rcParams['font.family'] = 'Nimbus Roman'
    plt.figure(figsize=(10,2))
    plt.plot(TDATE, TDATA, label='窗值', color = 'blue')
    plt.plot(DDATE, DDATA, marker = 'o', linestyle = '-', color = 'r', label = '日值', markersize=2)
    plt.title('Apprent displacement')
    plt.ylabel('Disp/um', fontsize=12, font='Nimbus Roman')
    plt.title('Time Series Plot')
    sdate, edate, xtick4Index, xtick4Label = sdate2edate4timeserial(TDATE)
    plt.xticks(xtick4Index, xtick4Label)  # 设置X轴刻度位置和标签
    plt.legend()
    plt.ylim(0, 0.3)
    plt.savefig(ouRoot + 'output.jpg', dpi=300, bbox_inches='tight')
    # 显示图表（可选）
    print(ouRoot + 'output.jpg')
    plt.show()  
# 示例数据
def paintdisp4dtvalue(TDATE, TDATA):
    #TDATE = ['2023-01-01', '2023-01-02', '2023-01-03', '2023-01-04']
    #TDATA = [100, 200, 150, 250]

    TDATA = TDATA;
    # 绘制图表
    plt.figure(figsize=(10, 6))
    plt.plot(TDATE, TDATA, color='blue',label="位移")
    plt.xlabel('Date[yy-mm-dd]', fontsize=10)
    plt.ylabel('Disp/um', fontsize=12)
    plt.title('Time Series Plot')
    sdate, edate, xtick4Index, xtick4Label = sdate2edate4timeserial(TDATE)
    plt.xticks(xtick4Index, xtick4Label)  # 设置X轴刻度位置和标签

    # 保存图表
    plt.savefig('output.jpg', dpi=300, bbox_inches='tight')
    # 显示图表（可选）
    plt.show()

def sdate2edate4timeserial(T):
    row, col     = T.shape;
    sdate        = np.round(T[0]);
    edate        = np.round(T[row - 1]);
    xtick4Index  = np.arange(sdate, edate, 5);
    N            = len(xtick4Index);
    xtick4Label  = [""] * N
    for i in range(0, N):
        time_obj       = py_tool.python_datetime(xtick4Index[i]);
        xtick4Label[i] = py_tool.datestr(time_obj, '%y %m-%d')
    return sdate, edate, xtick4Index, xtick4Label
def sdate2edate4timeserial4any(T):
    row, col     = T.shape;
    sdate        = np.round(T[0]);
    edate        = np.round(T[row - 1]);
    #print(sdate)
    #print(edate)
    s_dt = float(sdate[0])
    e_dt = float(edate[0])
    time_obj        = py_tool.python_datetime(s_dt)
    [iy,im,id,iH,iM,iS] = py_tool.datevec(time_obj)
    #print(f"{iy}-{im}-{id}")
    time_obj        = py_tool.python_datetime(e_dt)
    [jy,jm,jd,jH,jM,jS] = py_tool.datevec(time_obj)
    #print(f"{jy}-{jm}-{jd}")
    sdate = []
    edate = []
    xtick4Index = []
    xtick4Label = []
    if jd > 0:
       jy = jy + 1
    if (jy - iy) >= 4:
       sdate, edate, xtick4Index, xtick4Label = calLabel4Year(iy, jy)
    else:
       sdate, edate, xtick4Index, xtick4Label = calLabel4Month(iy, im, id, jy, jm, jd)
    xtick4Month = calTick4Month(iy, jy)
    return sdate, edate, xtick4Index, xtick4Label, xtick4Month

def calTick4Month(iy, jy):
    xtick4year   = np.arange(iy, jy + 1, 1)
    N            = len(xtick4year)
    xtick4Index  = np.full(N * 12, np.nan)
    for i in range(0, N):
        cy = i + iy
        for j in range(0, 12):
            cm = j + 1
            xtick4Index[i * 12 + j] = py_tool.python_datenum(cy, cm, 1) # matlab date
    return xtick4Index
def calLabel4Year(iy, jy):
    print(iy)
    print(jy)
    xtick4year   = np.arange(iy, jy + 1, 1)
    N            = len(xtick4year)
    print(N)
    xtick4Label  = [""] * N
    xtick4Index= np.full(N, np.nan)
    ind        = 0
    sdate   = py_tool.python_datenum(iy, 1, 1) #matlab date
    edate   = py_tool.python_datenum(jy, 1, 1) #matlab date
    for i in range(0, N):
        cy = i + iy
        xtick4Index[i]   = py_tool.python_datenum(cy, 1, 1) # matlab date
        time_obj         = py_tool.python_datetime(xtick4Index[i]) #python date
        xtick4Label[i]   = py_tool.datestr(time_obj, '%Y')
        ind = ind + 1
    return sdate, edate, xtick4Index, xtick4Label

def calLabel4Month(iy, im, id, jy, jm, jd):
    #起止时间的设定
    sdate      = py_tool.python_datenum(iy, im, 1) #matlab date
    edate      = py_tool.python_datenum(jy, jm, 1) #matlab date
    if jd > 0:
       edate = py_tool.python_addtodate(edate, 1, 'months')

    #循环产生多少个标签
    cdate      = sdate
    N          = 0
    while cdate <= edate:
          N = N + 1
          cdate = py_tool.python_addtodate(cdate, 1, 'months')

    #循环产生标签
    xtick4Index = np.full(N, np.nan)
    xtick4Label  = [""] * N
    cdate        = sdate
    for i in range(0 , N):
        xtick4Index[i] = cdate
        time_obj       = py_tool.python_datetime(cdate)
        xtick4Label[i] = py_tool.datestr(time_obj, '%y-%m')
        cdate          = py_tool.python_addtodate(cdate, 1, 'months')

        
    return sdate, edate, xtick4Index, xtick4Label

def calYLim4Auto(T, alpha = 0.95):
    lo,up = py_tool.getULbynormal(T, alpha)
    return up, lo

def calYLim4Auto2MaxMin(T):
    max4u = np.nanmax(T)
    min4l = np.nanmin(T)
    dis   = max4u - min4l
    dis5  = dis / 5.0
    up    = max4u + dis5
    lo    = min4l - dis5
    return up, lo

if __name__ == "__main__":
    SDATE = 0
    EDATE = 3650
    EQRANK= [[]]
    calpaintEqLabelIndex(SDATE, EDATE, EQDATE, EQRANK)

    

        

    




